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HIGHLIGHTS 


►  A  reduced  order  model  is  proposed  for  solid-phase  diffusion. 

►  The  reduced  order  model  is  faster  than  solving  the  full  model  and  is  very  accurate. 

►  Impact  of  weight  function  on  solution  is  investigated. 

►  Step  response  method  is  proposed  for  low  frequencies. 

►  Complex  exponential  method  is  proposed  for  high  frequencies. 
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A  model  order  reduction  method  is  developed  and  applied  to  the  solid-phase  diffusion  problem  used  in 
physics-based  lithium  ion  cell  models.  The  reduced  order  model  is  in  the  form  of  a  state  space  model. 
Model  identification  is  performed  in  the  frequency-domain  using  the  vector  fitting  method.  The  method 
allows  the  user  to  control  the  order  of  the  model,  the  frequency  band  for  model  identification,  and 
optionally  a  weight  function  to  give  a  certain  frequency  band  more  weight.  The  model  can  be  used  for 
spherical  and  non-spherical  particles.  For  spherical  particles,  the  results  from  using  the  reduced  order 
model  are  compared  with  those  from  analytical  solutions,  and  excellent  agreement  is  achieved  using  3rd 
and  5th  order  models.  When  the  approach  is  applied  to  non-spherical  particles,  the  transfer  functions 
need  to  be  calculated  numerically.  Two  methods,  step  response  and  complex  exponential,  are  proposed 
to  calculate  the  required  transfer  function.  While  the  step  response  method  is  more  suitable  for  low 
frequencies,  the  exponential  method  is  more  accurate  for  high  frequencies. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  lithium  ion  battery  is  a  preferred  candidate  as  a  power 
source  for  hybrid  electric  vehicles  and  electric  vehicles  due  to  its 
outstanding  characteristics  such  as  high  energy  density,  high 
voltage,  low  self-discharge  rate,  and  good  stability  among  others. 
Physics-based  lithium  ion  battery  models  are  widely  used  to 
predict  the  electrochemical  behavior  of  lithium  ion  batteries  [1—3], 
The  implementation  of  such  a  model  requires  solving  a  diffusion 
problem  in  the  porous  electrodes.  The  diffusion  process  in  the 
electrodes  is  by  nature  3D  because  of  the  rather  convoluted  shape 
of  the  porous  structures.  Therefore,  it  is  computationally  expensive 
to  solve  the  physics-based  lithium  ion  battery  models  when  3D 
diffusion  equations  are  used  to  characterize  the  diffusion  process. 
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Due  to  this  difficulty,  the  porous  structure  is  assumed  to  consist  of 
disconnected  spherical  particles  so  that  the  diffusion  equation  can 
be  assumed  to  be  one-dimensional.  In  such  an  approach,  each 
particle  is  exposed  to  different  boundary  conditions,  and  therefore 
the  diffusion  equation  has  to  be  solved  for  each  particle.  Due  to  the 
large  number  of  particles  involved,  it  is  still  computationally 
expensive  to  solve  the  physics-based  lithium  ion  battery  models 
when  a  large  number  of  ID  diffusion  equations  are  solved 
numerically. 

Under  the  assumption  of  spherical  particles,  several  methods 
have  been  proposed  to  reduce  the  size  of  the  problem.  They  belong 
to  time-domain  approach  or  frequency-domain  approach.  One 
time-domain  approach  makes  the  assumption  that  the  concentra¬ 
tion  within  each  spherical  particle  can  be  approximated  with 
a  parabolic  profile  [4,5],  Another  time-domain  approach  obtains  an 
approximate  solution  by  truncating  the  analytical  solution  of  an 
infinite  series  and  then  adding  an  estimate  term  for  the  truncation 
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error  [6],  In  these  two  methods,  spherical  particles  are  assumed, 
and  the  methods  cannot  be  extended  to  non-spherical  particles. 
Time-domain  methods  also  cannot  control  frequency  performance 
easily,  while  frequency-domain  methods  can  achieve  this  readily.  In 
frequency-domain  approximation  methods,  a  state  space  model  is 
created  to  approximate  the  transfer  function  of  the  system.  In 
Ref.  [7],  high-order  poles  of  the  transfer  function  are  truncated  and 
the  lower-order  poles  are  grouped  together  and  approximated 
using  a  state  space  model.  The  proposed  5th  order  reduced  model 
gives  an  error  of  6.3%.  More  accurate  results  can  be  obtained  using 
the  method  discussed  in  this  paper.  In  Ref.  [8],  a  discrete-time  state 
space  model  is  derived  from  a  known  transfer  function.  And  then  it 
is  converted  into  a  continuous-time  state  space  model.  Using  the 
method  in  this  paper,  a  continuous-time  state  space  model  can  be 
generated  directly.  In  both  Refs.  [7,8],  analytical  transfer  functions 
are  assumed  to  be  known,  and  thus  both  methods  are  limited  to 
spherical  particles.  The  method  discussed  in  this  paper  can  be 
extended  to  particles  of  any  shape. 

In  Ref.  [9],  Hu  et  al.  proposed  a  continuous-time  state  space 
model,  which  is  created  directly  from  a  numerically  calculated 
transfer  function.  The  model  is  demonstrated  to  be  applicable  to 
particles  of  any  shape.  In  this  paper,  the  method  used  in  Ref.  [9]  is 
improved  and  refined.  It  was  found  that  the  step  response  method 
used  in  Ref.  [9]  to  calculate  the  transfer  function  is  very  accurate 
only  for  relatively  low  frequencies.  This  is  because  magnitudes  of 
high  frequency  components  are  much  smaller  than  those  of  lower 
frequency  components  and  therefore  the  step  response  method 
gives  large  numerical  error  for  high  frequency  components.  In  this 
paper,  a  complex  exponential  method  is  proposed  to  calculate  the 
transfer  function  for  high  frequencies.  It  will  be  shown  in  this  paper 
that  the  complex  exponential  method  gives  highly  accurate  results 
for  high  frequencies.  But  this  method  is  hard  to  converge  for  lower 
frequencies.  When  the  complex  exponential  method  is  used 
together  with  the  step  response  method,  accurate  transfer  function 
can  be  obtained  for  the  entire  interested  frequency  range.  Model 
identification,  after  the  transfer  function  is  determined,  is  done  in 
the  frequency-domain  using  the  vector  fitting  (VF)  method  [10]. 
The  VF  method  is  a  widely  used  method  for  model  identification 
[11—16],  and  it  is  the  method  used  in  Ref.  [9],  In  this  paper,  a  weight 
function  is  added  during  the  fitting  process  and  its  impact  on  the 
model  performance  is  investigated. 

2.  Model  development  when  analytical  transfer  functions  are 
available 

In  using  the  frequency-domain  approach,  the  problem  is  treated 
like  a  system.  A  system  is  an  entity  that  processes  a  set  of  input 
signals  (or  simply  called  inputs)  and  yields  another  set  of  output 
signals  (or  simply  called  outputs).  In  such  a  system  view,  only  the 
input/output  relationship  of  the  system  is  of  interest  to  the  user  and 
the  inner  structure  of  the  system  is  not  important.  In  the  solid-phase 
diffusion  problem,  the  molar  flux  at  the  particle  surface  as  a  function 
of  time  is  used  as  the  input  and  the  difference  between  the  surface 
concentration  and  the  average  concentration  at  the  particle  surface 
as  a  function  of  time  is  used  as  the  output.  We  are  interested  in  the 
relationship  between  the  surface  molar  flux  and  the  difference 
between  the  surface  concentration  and  the  average  concentration. 
Concentration  distribution  inside  a  particle  is  not  of  great  interest. 
Solid-phase  diffusion  problem,  as  described  above,  is  a  system.  More 
importantly,  it  is  not  only  a  system,  but  it  is  also  a  linear  and  time- 
invariant  (LTI)  system.  Any  state  space  model  is  also  an  LTI  system. 
One  important  feature  about  LTI  systems  is  that  if  two  LTI  systems 
have  the  same  transfer  function,  the  two  systems  behave  identically 
in  that  the  outputs  of  the  two  systems  are  the  same  provided  that  the 
inputs  to  the  two  systems  are  the  same.  This  feature  allows  us  to  use 


the  state  space  model  to  simulate  the  solid-phase  diffusion  problem 
provided  that  its  transfer  function  is  curve-fitted  to  that  of  the  solid- 
phase  diffusion  problem.  Since  the  size  of  the  state  space  model  is 
small,  it  runs  very  fast  compared  with  a  full  model  which  solves  the 
complete  diffusion  equation  directly  using  numerical  methods. 

To  demonstrate  the  state  space  approach,  we  first  model 
a  spherical  particle,  for  which  an  analytical  transfer  function  is 
available.  The  diffusion  of  lithium  ions  in  a  spherical  particle 
follows  the  Fick’s  law  and  is  described  by  the  following  partial 
differential  equations: 


c(r,  t  =  0)  =  0 


(2) 


(3) 


=m  (4) 

drlr=R 

where  c  is  the  concentration  of  lithium  in  the  solid  particle;  D  is  the 
solid-phase  diffusion  coefficient  for  lithium  in  the  particle,  assumed 
to  be  constant;  j(t)  is  the  boundary  molar  flux,  assumed  to  be 
a  function  of  time;  r  is  the  radial  coordinate;  R  is  the  radius  of  the 
spherical  particle.  Note  that  we  specified  the  initial  condition  to  be 
zero  in  Eq.  (2)  because  such  an  initial  condition  is  convenient  for 
Laplace  transform.  The  real  initial  value  can  simply  be  added  back 
to  the  final  solution.  We  are  interested  in  the  relationship  between 
j(t)  and  surface  concentration,  cs  =  c(R,  t).  From  Eqs.  (1)— (4), 
Jacobsen  and  West  [17]  found  the  transfer  function: 

Cs(s)  R  [  tanh(0)  ] 

M  -  D  [tanh<0)  -  0\  l  J 

with  =  Ry/sjD.  We  define  the  difference  between  the  surface 
concentration  and  the  average  concentration  by  subtracting  the 
average  response  Cave(s)  =  -3J(s)/(Rs)  from  the  surface  concen¬ 
tration.  Cave  (s)  is  obtained  by  taking  the  Laplace  transform  of  the 
mass  conservation  equation,  (4/3)7rR3(dcave(t)/dt)  =  j(t) -4nR2. 
Define  ACs(s)  =  Cs(s)  -  Cave(s),  and  the  transfer  function  relating 
ACs(s)  to  J(s)  becomes  the  following: 


ACs(s)  R  l"(^2  +  3)tanh(/J)  -  3/Jj 

Ks)  ~  D[  /?2(tanh(/?)  -  0)  J 


(6) 


An  analytical  solution  of  Eqs.  (1 )— (4)  under  a  unit  step  input  in 
non-dimensional  form  is  shown  to  be  an  infinite  series  in  Ref.  [6], 
After  converting  it  into  dimensional  form,  the  surface  concentra¬ 
tion  relative  to  the  average  concentration  under  a  unit  step  input  of 
1  mol  m-2  s-1  is  the  following, 


ACS, Step (t)  =  E  Qn(f)  £  0)  (7) 

where 

Qn(t)  =  “s|[1“eXP(“A"^)]  W 

A„-tan(An)  =  0  n  =  l,2,...  (9) 

The  analytical  step  response  shown  above  will  be  used  to  vali¬ 
date  models  in  the  time-domain.  Note  that  A Cs(s)  in  Eq.  (6)  is  the 
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Laplace  transform  of  A cs(t).  AcSiStep(t)  in  Eq.  (7)  is  the  step  response 
of  surface  concentration  relative  to  the  average  concentration  in  the 
time-domain. 

A  general  state  space  model  is  typically  written  as  follows: 


x  =  Ax  +  Bu 
y  =  Cx  +  Du 


(10) 


x  is  the  state  vector;  y  is  the  output  vector;  u  is  the  input  vector;  A,  B, 
C,  D  are  constant  coefficient  matrices  of  proper  sizes.  Since  Eq.  (10) 
is  linear  with  constant  coefficients,  the  system  is  an  LTI  system.  For 
the  state  space  model  used  to  simulate  the  solid-phase  diffusion 
problem,  x  has  no  physical  meaning;  y  is  a  scalar  representing 
surface  concentration  relative  to  the  average  concentration,  Acs(t); 
u  is  a  scalar  representing  the  transient  wall  flux,  j(t).  Matrix  D 
becomes  zero  for  the  diffusion  problem  since  a  step  change  in  j(t) 
does  not  cause  a  step  increase  in  surface  concentration.  Taking  the 
Laplace  transform  of  Eq.  (10)  gives  the  transfer  function  of  the  state 
space  model, 


Y(s)/U(s)  =  CisI-Ay'B 

Its  equivalent  rational  form  is  the  following: 


HC0(c,a,s)  =  Y(s)/U(s )  = 


(11) 


(12) 


function  already  gives  very  negligible  steady  state  error.  So, 
a  weight  function  is  hardly  needed  to  minimize  the  steady  state 
error  using  the  VF  method.  Weight  functions  will  be  discussed  in 
details  in  Section  3. 

In  using  the  state  space  model  to  model  the  diffusion  process, 
the  goal  is  to  identify  Eq.  (12)  with  the  system  transfer  function,  Eq. 
(6),  for  s  =jw.  If  we  denote  the  sampled  system  transfer  function  of 
Eq.  (6)  to  be  H(Sj  =  jw{)  i  =  1  ,...,M,  with  M  >  N,  the  problem 
amounts  to  generating  a  state  space  model  of  Eq.  (10)  such  that  the 


E(c,  a,s)  =  Hca(  c,  a,s)  -  H(s)  (17) 

evaluated  at  the  discrete  frequency  points,  is  minimized  in  some 
appropriate  norm.  A  least-squares  solution  can  be  obtained  by 
minimizing 


f(c ,  a)  =  aJ  E  IIWcq(  c ’  a , s,-)  —  H(s,-)||2  (18) 

Minimization  of  Eq.  (18)  is  a  non-linear  problem.  The  VF  method 
is  used  in  this  paper,  which  is  a  robust  iterative  method  to  solve  the 
minimization  problem.  A  brief  explanation  of  the  VF  method 
applied  to  the  diffusion  problem  is  shown  in  Appendix  B.  Details 
about  the  VF  method  can  be  found  in  Ref.  [10]. 


where  c  =  [ci,  ...,qv],  a  —  [a\ a^]  are  the  residuals  and  poles 
of  the  transfer  function,  and  N  is  the  order  of  the  state  space  model. 

From  the  final-value  theorem,  we  could  put  a  constraint  on  Eq. 
(12),  namely  lim  Eq.(12)  =  lim  Eq.(6)  (see  Appendix  A),  to 
enforce  that  Eq)X*fc)  has  no  steacfy  state  error, 


E-- 

ai 


R 

5D 


(13) 


An  equivalent  form  of  Eq.  (12)  with  the  constraint  of  Eq.  (13)  is 
the  following, 


Y{s)/U(s)  =  Z  +  Es!!~.  (14) 

subject  to, 

Z+^c,  =  0  (15) 


where  q  and  a,  are  new  model  constants  similar  to  q  and  a,  that 
need  to  be  identified  if  the  constrained  version  of  the  transfer 
function  of  Eq.  (14)  is  used.  An  advantage  of  Eq.  (14)  is  that  it  gives 
an  explicit  expression  for  the  steady  state  value,  namely  Z.  The 
constraint  of  Eq.  (15)  is  to  satisfy  initial  value  of  zero  under  unit  step 
input.  This  is  derived  in  a  similar  fashion  as  shown  in  Appendix  A 
using  the  initial-value  theorem. 

Eq.  (12)  without  the  constraint  of  Eq.  (13)  is  used  in  this  paper 
mainly  because  the  VF  method  used  in  this  paper  requires  the  form 
of  Eq.  (12).  However,  the  method  does  allow  one  to  minimize  the 
steady  state  error  by  using  a  different  approach  rather  than  using 
the  constraint  of  Eq.  (13).  More  specifically,  when  a  weight  function 
with  a  heavy  weight  on  frequency  of  zero  is  used,  the  steady  state 
error  will  be  minimized.  Note  that  the  VF  method  without  a  weight 


2.1.  A  spherical  particle  test 

Property  values  used  in  the  test  case,  defined  in  Table  1,  are 
typical  of  solid  state  diffusion  in  electrochemical  cells.  The  char¬ 
acteristic  time  t~(R2/D)  =  5000  s  indicates  that  it  can  take  over  an 
hour  for  solid-phase  concentration  gradients  to  relax.  High  power 
HEV  batteries,  however,  may  become  solid  state  transport  limited 
in  as  little  as  5  s  [18],  To  capture  these  disparate  timescales,  solution 
up  to  10  Hz  is  sampled  for  model  identification.  Such  a  frequency  is 
referred  to  as  cut-off  frequency  for  model  identification  in  this 
paper. 

Fig.  1  compares  the  frequency  response  of  the  state  space 
models  with  the  analytical  frequency  response,  Eq.  (6).  The  5th 
order  state  space  model  gives  better  accuracy  than  the  3rd  order 
state  space  model.  Results  in  Fig.  1  do  not  appear  very  accurate 
especially  for  high  frequency  components.  This  is  because  the  VF 
method  by  default  uses  the  2-norm  to  minimize  the  error  in  the 
frequency-domain  without  a  weight  function,  as  shown  in  Eq.  (18). 
Since  high  frequency  components  have  much  less  magnitudes,  as 
shown  in  Fig.  1,  VF  sacrifices  their  accuracy.  If  one  checks  the  error 
quantitatively,  the  VF  error  in  2-norm  is  actually  only  3.2%  and 
0.67%  for  the  3rd  order  and  5th  order  model,  respectively.  The  large 
apparent  error  in  Fig.  1  is  due  to  logarithmic  scale  used  in  the  plot. 
Fig.  2  compares  the  solution  in  time-domain  under  unit  step  input. 
The  “analytical”  solution  in  the  time-domain  is  approximated  by 
taking  the  first  5000  terms  in  Eq.  (7).  Both  models  give  excellent 
results  compared  with  the  analytical  solution  as  shown  in  Fig.  2. 
The  small  difference  between  the  3rd  order  state  space  model  and 
the  analytical  solution  is  only  visible  in  the  logarithmic  time  scale 
plot,  Fig.  2b.  Most  signals,  including  the  unit  step,  have  large 
portion  of  low  frequency  components  and  small  portion  of  high 


Spherical  particle  model  parameters. 
Parameter  value 
Diffusion  coefficient,  D  (m2  s  ') 
Particle  radius,  R  (m) 


Value 

2.0  x  10-16 
1.0  x  10~6 
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Fig.  1.  Frequency  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with 
those  from  the  3rd  order  and  the  5th  order  space  models. 


Time  (sec) 


Time  (sec) 

Fig.  2.  Unit  step  responses  from  the  analytical  solution  and  two  state  space  models,  a) 
Regular  time  scale,  b)  Log  time  scale. 


frequency  components.  So,  the  2-norm  used  by  VF  without 
a  weight  function  makes  time-domain  results  accurate  for  most 
signals.  This  explains  very  accurate  time-domain  results  by  just 
a  3rd  order  state  space  model.  Also  note  that  results  in  Fig.  2  show 
negligible  steady  state  error  without  using  any  weight  function. 

If  one  has  a  signal  with  only  high  frequency  components,  then 
even  the  5th  order  model  will  not  give  accurate  results.  To 
demonstrate  this,  a  10th  order  model  is  generated,  which  gives 
much  better  accuracy  in  high  frequency  band  compared  with  the 
5th  order  model.  The  frequency  response  of  the  10th  model  is 
shown  in  Fig.  3.  A  pure  sinusoidal  signal  of  1  Hz  is  sent  to  the  5th 
order  state  space  model  and  the  10th  order  state  space  model.  Their 
responses  are  compared  in  Fig.  4.  It  is  clear  the  5th  order  state  space 
model  deviates  from  the  solution  of  the  10th  order  state  space 
model.  But  practically,  when  the  reduced  model  is  integrated  into 
the  physics-based  cell  models,  the  reduced  model  is  unlikely  to 
experience  such  a  high  single  frequency  signal.  Even  if  such  a  single 
frequency  signal  is  used  at  the  boundary  of  the  physics-based 
model,  the  particle  wall  flux  will  not  have  such  a  single  frequency 
since  the  physics-based  model  is  non-linear.  However,  if  significant 
high  frequency  components  are  present,  VF  can  generate  higher 
order  models  to  give  accurate  results.  In  such  a  case,  the  full  model 
also  needs  more  grid  resolution  for  accuracy.  So,  it  is  still  beneficial 
to  use  the  reduced  model. 


3.  Vector  fitting  with  weight  functions 

A  weight  function  can  be  added  when  performing  VF  to  give 
more  weight  to  certain  frequency  band.  Eq.  (18)  with  a  weight 
function  becomes: 


f(c,  a)  =  ^f>(s;H|Hca(c,  a,Sj)  -H(S() | |j  (19) 

A  weight  function  giving  high  frequency  components  more 
weight  is  used  to  demonstrate  the  effect  of  the  weight  function  on 
the  frequency  response.  Since  the  magnitude  of  the  system 
response  decreases  exponentially,  the  weight  function  used  has  an 
exponential  form  of  1.07'  to  bring  all  frequency  components  to 
approximately  the  same  magnitude.  The  number  1.07  was  chosen 
by  trial  and  error  to  give  the  results  shown  in  Fig.  5.  It  is  clear  from 
Fig.  5  that  the  frequency  response  from  the  model  with  a  weight 
function  favoring  high  frequencies  appears  to  be  more  accurate. 
However,  the  time-domain  performance  of  the  model  under  a  unit 
step  input  is  worse  as  shown  in  Fig.  6.  This  is  because  while  favoring 
high  frequency  its  low  frequency  accuracy  is  compromised.  A  unit 
step  signal,  like  many  other  signals,  has  more  components  in  low 
frequency  than  in  high  frequency.  Therefore,  the  time-domain 
performance  of  a  model  favoring  high  frequencies  is  not  as  accu¬ 
rate  under  a  unit  step  input.  Note  that  the  weighted  5th  order  state 
space  model  has  a  noticeable  steady  state  error.  Using  the  final- 
value  theorem,  the  same  conclusion  can  be  drawn  from  its 
frequency  response  as  frequency  approaches  zero.  A  way  to 
improve  the  steady  state  solution  is  to  put  more  weight  on  very  low 
frequency  components.  But  if  we  put  weight  on  both  low  and  high 
frequencies,  the  midfrequency  components  will  have  poor  perfor¬ 
mance.  After  all,  we  only  have  10  degrees  of  freedom  at  our  disposal 
to  optimize  Eq.  (19)  using  a  5th  order  state  space  model. 

Ultimately  the  state  space  model  will  be  integrated  into  the 
physics-based  cell  models  and  these  models  are  highly  non-linear 
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Fig.  3.  Frequency  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with  pig.  5.  Frequency  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with 
those  from  the  10th  order  state  space  model.  those  from  two  5th  order  state  space  models  with  and  without  weight  function. 


overall.  As  shown  above,  adding  weight  function  to  make  frequency 
response  look  more  accurate  might  not  always  be  desirable.  To 
make  sure  that  the  state  space  model  gives  accurate  results  for 
enough  frequency  range,  one  could  run  the  full  physics-based 
model  using  signals  containing  reasonably  high  frequency 
components  and  record  the  j(t)  at  the  particle  surface.  An  FFT  ofj(t) 
tells  what  frequency  contents  are  present  at  the  particle  surface. 
Based  on  that  information,  a  judicious  decision  can  be  made  as  for 
what  cut-off  frequency  to  be  used  and  what  weight  function,  if  any, 
to  be  used  when  generating  the  state  space  model.  The  benefit  of 
using  the  VF  method  in  identifying  state  space  models  is  that  it  can 
accommodate  all  these  requests  easily  and  robustly.  More  details 


on  how  to  integrate  the  reduced  order  model  into  the  physics- 
based  cell  models  are  shown  in  Ref.  [9], 

4.  Model  development  when  analytical  transfer  functions  are 
not  available 

We  showed  how  to  create  state  space  models  for  particle 
diffusion  problems  when  the  transfer  function  is  known.  For 
spherical  particles,  the  transfer  function  of  the  system  has  an 
analytical  form  of  Eq.  (6).  For  non-spherical  particles,  we  need  to 
approximate  system  transfer  functions  numerically  before  VF  can 
be  used  to  identify  state  space  models.  In  this  section,  we  demon¬ 
strate  how  to  obtain  transfer  functions  accurately  using  numerical 


Fig.  4.  Responses  of  5th  order  and  10th  order  state  space  models  under  1  Hz  sinusoidal  Fig.  6.  Step  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with  those  from 
signal.  Solid  line  is  for  10th  order;  dash  line  is  for  5th  order.  two  5th  state  space  models  with  and  without  weight  function. 
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methods.  Spherical  particles  will  be  used  since  analytical  solution 
exists  to  help  validate  the  methodology.  So,  the  analytical  transfer 
function  is  not  used  in  the  model  identification  process.  The 
analytical  transfer  function  is  only  used  to  validate  results. 

4.1.  Step  response  method  to  obtain  frequency  response 

One  way  to  calculate  the  system  transfer  function  is  to  calculate  the 
Laplace  transform  of  the  impulse  response  of  the  system.  As  a  matter 
of  fact,  Laplace  transform  of  the  impulse  response  of  a  system  is  one 
way  to  define  transfer  function.  But  since  it  is  difficult  to  obtain  the 
impulse  response  accurately,  we  obtain  the  step  response  first  and 
then  calculate  the  time  derivative  of  the  step  response  to  obtain  the 
impulse  response.  Such  a  procedure  can  be  used  for  non-spherical 
particles  since  step  responses  can  be  obtained  by  numerical 
methods  regardless  of  the  shape  of  the  particle. 

The  response  of  a  spherical  particle  under  a  unit  step  input  of 
1  mol  m-2  s-1  is  Eq.  (7).  The  time  derivative  of  the  unit  step 
response  becomes  the  unit  impulse  response,  which  is  the 
following, 

where  A„  is  defined  in  Eq.  (9).  If  the  unit  step  response  is  obtained 
by  numerical  methods,  the  unit  impulse  response  can  be  calculated 
by  taking  numerical  derivative. 

Fourier  transfer  of  Eq.  (20)  becomes  the  frequency  response  of 
the  system,  namely  Eq.  (6)  with  s  =  jui.  However,  to  make  the 
process  general,  we  are  not  assuming  we  could  obtain  analytical 
expression  for  the  impulse  response.  So,  sampled  version  of  Eq.  (20) 
will  be  used.  FFT  of  the  sampled  impulse  response,  after  proper 
scaling,  becomes  the  numerically  calculated  frequency  response  at 
discrete  frequency  values.  And  frequency  response  at  discrete 
frequency  values  is  what  we  need  to  perform  VF  to  identify  the 
state  space  model.  Fig.  7  shows  three  numerically  obtained 
frequency  responses  by  using  sampling  frequencies  of  5  Hz,  50  Hz, 
and  100  Hz,  respectively.  As  sampling  frequency  becomes  higher, 
the  numerically  calculated  frequency  response  becomes  more 
accurate.  Note  that  we  obtained  the  sampled  impulse  response 
from  the  analytical  solution  of  Eq.  (20),  approximated  by  taking  the 
first  5000  terms.  So,  there  is  a  minimum  amount  of  numerical  error 
in  this  process.  When  an  impulse  response  is  calculated  numeri¬ 
cally  from  a  step  response,  which  in  turn  is  from  numerical 
methods,  it  is  difficult  to  obtain  accurate  results  for  high  frequen¬ 
cies.  Therefore,  the  step  response  method  can  only  give  accurate 
results  to  relatively  low  frequencies. 

4.2.  Complex  exponential  method  to  obtain  frequency  response 

The  step  response  method  can  be  used  to  obtain  accurate 
frequency  response  for  low  frequencies.  But  it  could  involve  large 
numerical  error  for  high  frequencies.  The  complex  exponential 
method  proposed  in  this  section  can  be  used  to  obtain  accurate 
frequency  response  at  discrete  frequency  values  for  high  frequen¬ 
cies.  Such  a  method  is  widely  used  in  high  frequency  electromag¬ 
netics  to  calculate  the  solution  in  the  frequency-domain  [19],  The 
complex  exponential  method  relies  on  the  fact  when  a  complex 
exponential  function  exp(/w0t)  is  used  as  an  input  to  a  linear 
system,  the  complex  amplitude  of  the  output  equals  the  frequency 
response  of  the  system  at  w0.  This  can  be  derived  from  the  fact  that 
the  Fourier  transfer  of  the  complex  exponential  is  2to5(w  -  w0). 
where  S  is  a  Dirac  delta  function.  This  feature  allows  us  to  calculate 
the  frequency  response  at  wo  by  calculating  the  complex  amplitude 
of  the  output  for  an  input  of  exp(/wot). 


1.0E-04  1.0E-03  1.0E-02  1.0E-01  1.0E+00  1.0E+01 

f  (Hz) 

Fig.  7.  Frequency  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with 
those  from  the  step  response  method. 


Let  j(t)  =  exp(jojt)  and  c(r,  t)  =  C(r)exp(j<,jf),  where  C  is 
a  complex  concentration  distribution.  Substituting  the  above  into 
Eq.  (1)  and  separating  real  and  complex  parts  gives  the  following 
governing  equation  for  complex  concentration  distribution,  C, 


m  m 


(21) 


where  Cr  and  Q  are  the  real  and  imaginary  part  of  C.  Note  that  Eq. 
(21 )  has  no  time  dependency.  The  frequency  response  we  need  is 
for  the  surface  concentration  relative  to  average  concentration.  The 
complex  average  concentration  from  complex  exponential  input 
can  be  shown  to  be  3j/(Ro>).  After  solving  Eq.  (21 )  for  C,  evaluating  C 
at  the  surface  minus  complex  average  concentration  gives  the 
desired  frequency  response  at  frequency  of  w.  Solving  Eq.  (21) 
multiple  times  for  discrete  values  of  w  gives  the  frequency 
response  at  those  discrete  frequency  values. 

Eq.  (21 )  is  a  pair  of  coupled  Poisson  equations,  which  can  be  solved 
easily.  FLUENT,  a  computational  fluid  dynamics  (CFD)  code  by  ANSYS 
Inc.  [20],  is  used  to  solve  Eq.  (21)  in  this  paper.  Fig.  8  shows  the 
frequency  response  at  discrete  frequency  values  calculated  using  the 
complex  exponential  method  compared  with  the  analytical  frequency 
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Fig.  8.  Frequency  responses  from  the  analytical  solution  of  Eq.  (6)  compared  with 
those  from  the  complex  exponential  method. 

response.  Fig.  8  shows  that  this  method  gives  highly  accurate  results 
for  frequency  as  high  as  10  Hz,  while  the  step  response  method,  when 
sampling  even  at  100  Hz,  is  less  accurate  at  high  frequencies  as  shown 
in  Fig.  7.  At  low  frequencies,  both  methods  give  excellent  results.  Note 
that  the  complex  exponential  method  solves  a  steady  state  problem 
rather  than  a  transient  problem.  Therefore,  in  using  the  complex 
exponential  method,  we  do  not  have  to  concern  about  time.  We  also 
do  not  need  to  worry  about  sampling  in  using  the  complex  expo¬ 
nential  method.  But  in  using  the  complex  exponential  method, 
multiple  steady  state  simulations  are  needed  as  opposed  to  solving 
just  one  transient  simulation  in  the  step  response  method.  When 
using  the  complex  exponential  method,  we  noticed  that  for  very  low 
frequencies  it  is  hard  to  converge  Eq.  (21 ).  So,  when  accurate  high 
frequency  solution  is  not  needed,  the  step  response  method  is  suffi¬ 
cient.  Otherwise,  the  complex  exponential  method  could  be  used  to 
calculate  the  frequency  response  at  high  frequencies. 

An  interesting  observation  was  obtained  when  solving  Eq.  (21). 
Fig.  9  shows  the  contour  of  magnitude  of  complex  concentration  at 
a  frequency  of  0.1  Hz.  Note  that  even  at  such  a  seemingly  low 
frequency  of  0.1  Hz  concentration  is  highly  concentrated  near  the 
surface.  This  suggests  that  the  curvature  of  the  sphere  has  negli¬ 
gible  impact  on  the  solution  at  this  frequency.  To  verify  this,  we 
utilize  the  analytical  solution  that  exists  when  the  radius  is  infinity, 
namely  a  wall  instead  of  a  sphere.  The  analytical  frequency 
response  can  be  shown  to  be  the  following  (see  Appendix  C), 


H(j(0)  = 


(22) 


The  above  frequency  response  is  compared  with  that  of  a  spherical 
particle  with  a  radius  of  1  pm  in  Fig.  10.  It  is  clear  for  Fig.  10  that  the 
behavior  of  the  spherical  particle  starts  to  differ  from  a  wall  at 
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Fig.  9.  Contour  of  magnitude  of  complex  concentration  C  (mol  m  3)  at  a  frequency  of 
0.1  Hz. 

a  frequency  of  approximately  0.1  Hz.  From  this  result,  it  is  clear  that 
0.1  Hz  for  this  particle  is  really  high. 

This  rather  interesting  observation  can  be  further  explained  by 
penetration  length,  defined  by  the  location  where  the  magnitude 
has  decreased  to  1%  of  the  wall  value.  For  the  wall,  the  penetration 
length  has  the  following  analytical  expression  (see  Appendix  C), 


f  (Hz) 


f  (Hz) 


Fig.  10.  Analytical  frequency  responses  of  Eq.  (6)  for  a  sphere  compared  with  those 
from  Eq.  (C.9)  for  a  wall. 
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Fig.  11.  Magnitude  of  complex  concentration  C  as  a  function  of  r  for  the  spherical 
particle  at  a  frequency  of  0.1  Hz  under  wall  flux  with  unit  magnitude. 

At  a  frequency  of  0.1  Hz,  the  penetration  length  is  approximately 
0.116  pm.  For  the  spherical  particle,  the  penetration  length  can  be 
found  by  plotting  the  concentration  as  a  function  of  r,  which  is  shown 
in  Fig.  11.  From  the  plot,  the  penetration  length  can  be  found  to  be 
0.121  pm.  This  is  within  5%  of  the  value  for  the  wall  confirming  that  the 
particle  with  a  radius  of  1  pm  behaves  like  a  wall  at  a  frequency  of 
0.1  Hz!  Such  information  helps  designers  to  determine  what  the 
optimum  particle  size  is  for  the  diffusion  process.  If  the  wall  flux  j(t) 
has  a  large  portion  of  frequency  components  close  to  0.1  Hz,  it  is 
clearly  advantageous  to  use  smaller  particles  to  increase  the  particle 
specific  interfacial  area  to  enhance  diffusion.  Note  that  the  above 
conclusion  is  based  on  the  diffusion  coefficient  used  in  Table  1. 

Once  the  reduced  order  model  is  identified  for  the  diffusion 
process,  the  model  can  be  integrated  into  the  physics-based  New¬ 
man  electrochemistry  pseudo-2D  model.  Such  an  integration 
progress  and  its  validation  against  solving  the  diffusion  equations 
numerically  are  discussed  in  Ref.  [9], 

5.  Conclusion 

A  model  order  reduction  method  is  developed  and  applied  to 
the  solid-phase  diffusion  problem  used  in  physics-based  lithium 
ion  cell  models.  Model  identification  is  performed  in  the  frequency- 
domain  using  the  vector  fitting  method.  The  method  allows  the 
user  to  control  the  order  of  the  model,  the  frequency  band  for 
model  identification,  and  optionally  a  weight  function  to  give 
certain  frequency  band  more  weight.  The  model  shows  excellent 
accuracy  with  a  5th  order  model.  The  3rd  order  model  shows 
a  small  deviation  from  analytical  solution  when  logarithmic  time 
scale  is  used.  Since  the  reduced  model  solves  only  a  few  equations, 
it  runs  much  faster  than  the  model  for  the  full  diffusion  equation. 

The  method  can  be  used  for  non-spherical  particles.  When 
doing  so,  system  responses  need  to  be  calculated  numerically.  Two 
methods  are  demonstrated.  The  step  response  method  is  easy  to 
use  and  gives  accurate  results  for  low  frequencies.  The  complex 
exponential  method  could  be  used  to  obtain  highly  accurate  results 
at  high  frequencies  but  it  is  hard  to  converge  at  low  frequencies. 

List  of  symbols 


c  concentration,  mol  m  3 

Cave  average  particle  concentration,  mol  m  3 

cs  surface  concentration,  mol  m~3 

C  complex  concentration 

Cave  Laplace  transform  of  average  concentration 


Cs  Laplace  transform  of  surface  concentration 
D  diffusivity  of  lithium  in  the  solid  particles,  m2  s-1 

R  radius  of  the  spherical  particle,  m 

j  wall  flux  of  lithium  ions,  mol  m  2  s-1 

J  Laplace  transform  of  wall  flux  of  lithium  ions 

w  weight  function  for  vector  fitting 

Acs  surface  concentration  relative  to  average  concentration, 
mol  m~3 

ACs.step  step  response  of  surface  concentration  relative  to  average 
concentration,  mol  m~3 

ACS  Laplace  transform  of  surface  concentration  relative  to 
average  concentration,  Cs  -  Cave 
3  penetration  length,  m 

An  eigenvalues  of  the  diffusion  problem 


Abbreviations 

FFT  Fast  Fourier  Transform 

LTI  Linear  and  Time-Invariant 

VF  Vector  Fitting 


Appendix  A.  A  constraint  on  numerical  transfer  function  to 
enforce  zero  steady  state  error  under  step  input  from  the 
final-value  theorem 


The  final-value  theorem  says  the  following.  If  A Cs(s)  is  the 
Laplace  transform  of  Acs(t),  then 

hm  Acs(t)  =  lim  sACs(s)  (A.1) 

From  Eq.  (6),  we  obtain  ACs(s)  =  H(s)J(s).  If  J(s)  is  the  Laplace 
transform  of  a  unit  step  function,  then, 

ACs(s)  =  (A.2) 

Substituting  Eqns.  (A.2)  into  (A.1)  gives  the  following, 


hm  Acs(t)  =  lim  H(s)  (A.3) 

Therefore,  if  we  want  to  enforce  zero  steady  state  error  under  step 
input,  we  need  the  following  constraint  on  the  numerical  transfer 
function,  namely, 


=  lim  H(s ) 


(A.4) 


Appendix  B.  Model  identification  process  for  the  diffusion 
problem  using  vector  fitting 

In  using  the  VF  method  for  state  space  model  identification,  the 
transfer  function  of  the  state  space  model  is  identified  with  the 
transfer  function  of  the  diffusion  problem.  The  transfer  function  of 
the  state  space  model  in  Eq.  (10)  can  be  shown  as  following: 

m  -  (b.i) 

The  fitting  process  can  be  viewed  as  using  rational  basis  func¬ 
tions  l/(s  —  a,)  for  curve-fitting  to  find  q.  If  a,-  were  known,  this 
would  become  a  linear  least-squares  problem  to  determine  q.  Since 
a,  are  not  known,  an  iterative  scheme  is  needed  to  update  a,  at  each 
iteration.  Once  the  poles  and  residuals  are  obtained,  it  is  a  simple 
realization  problem  to  obtain  the  corresponding  state  space  model 
of  Eq.  (10).  Such  a  fitting  approach  in  the  frequency-domain  is 
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generally  referred  to  as  rational  fitting.  VF  is  a  robust  algorithm  of 
rational  fitting  first  proposed  in  Ref.  [10]. 

The  steps  to  identify  the  state  space  model  using  VF  are 
summarized  as  follows: 

Step  1.  Obtain  the  transfer  function  of  the  diffusion  problem 
numerically  or  analytically. 

Step  2.  Perform  vector  fitting  to  obtain  the  poles  and  residuals  of 
the  transfer  function  of  the  state  space  model. 

Step  3.  Construct  the  state  space  model  from  the  known  transfer 
function,  which  is  obtained  from  Step  2. 


Appendix  C.  Analytical  solution  for  a  wall  diffusion  problem 
under  complex  exponential  flux  boundary  condition 

The  governing  equations  for  the  diffusion  problem  in  an  infinite 
wall  are 


dc 

dt 


(C.1) 


c(z,  t  =  0)  =  0 


(C.2) 


'm-m 


(C.3) 


where  z  is  the  coordinate  going  into  the  wall.  The  extra  minus  sign 
in  Eq.  (C.3)  compared  with  Eq.  (4)  is  due  to  the  convention  that 
positive  j(  t)  points  outside  of  the  domain. 

Let  j(t)  =  exp(j'«t)  and  c(z,t)  =  C(z)exp(/wt),  where  C  is 
complex  concentration  distribution.  Substituting  the  above  into  Eq. 
(C.1)  and  Eq.  (C.3)  gives  the  following, 


0 2_C 
dz2 


(C.4) 


-D^  =  -1  (C.5) 

9zlz=0 

A  trial  solution  of  the  form  C  =  exp(/7cz)  leads  to  the  condition: 
k2  =  -§  (C.6) 

This  leads  to  a  solution  of  C 

C  =  ct  exp  [  -  (1  +  j)  sJYdz]  (c?) 

The  second  exponential  term  drops  since  it  goes  to  infinity  as  z 
goes  to  infinity.  Boundary  condition  Eq.  (C.5)  can  be  used  to 
determine  the  constant  in  Eq.  (C.7).  The  final  expression  for  C  is  the 
following 


c-S5expb'Wiz] 

Evaluating  C  at  z  =  0  gives 


(C.8) 


C(  0)  = 


~t+j 

^/4np 


(C.9) 


Eq.  (C.9)  is  the  frequency  response  of  surface  concentration.  To 
obtain  the  frequency  response  of  the  surface  concentration  relative 
to  the  average  concentration,  we  need  to  subtract  the  average 
value.  But  since  the  wall  has  infinite  size,  the  average  value  is  zero. 

To  obtain  the  penetration  length,  defined  by  the  location  where 
the  magnitude  has  decreased  to  1%  of  the  wall  value,  we  let  the 
exponential  term  in  Eq.  (C.8)  delays  to  0.01, 


(C.10) 


(C.11) 
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